Introduction

Pour cette deuxième séance, nous allons apprendre à utiliser un planifateur d’itinéraire sur R afin de calculer des isochrones autour de chaque supermarché et ainsi pouvoir identifier les zones de chalandises.

1. Cartographier des isodistances

Une première information simple à obtenir mais très informative peut être obtenue en traçant des courbes d’isodistances autour de chaque supermarchés. On peut ainsi indentifier les zones non-deserviés.

télécharger le package en .zip en cliquant sur ce lien
https://cran.r-project.org/src/contrib/Archive/opentripplanner/opentripplanner_0.4.0.tar.gz

Et installer Java 8 sur votre ordinateur:
https://www.java.com/fr/download/
library(dplyr)
library(ggplot2)
library(sf)
library(readr)
library(readxl)
library(tmap)
library(tmaptools)
library(here)

otp_available <- requireNamespace("opentripplanner", quietly = TRUE)
if (!otp_available) {
  otp_pkg <- here("Cours", "opentripplanner_0.4.0.tar.gz")
  if (file.exists(otp_pkg)) {
    try(install.packages(otp_pkg, repos = NULL, type = "source"), silent = TRUE)
    otp_available <- requireNamespace("opentripplanner", quietly = TRUE)
  }
}
## Warning in install.packages(otp_pkg, repos = NULL, type = "source"): 'lib =
## "/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library"' is not
## writable
if (otp_available) {
  library(opentripplanner)
} else {
  message("Package 'opentripplanner' absent: les sections OTP/isochrones seront ignorees.")
}

# Chargement des données locales (mêmes sources que test.Rmd)
iris_shp <- st_read(
  dsn = here("Cours"),
  layer = "georef-france-iris-millesime"
)
## Reading layer `georef-france-iris-millesime' from data source 
##   `/Users/loic/Documents/UGA/donnee_entreprise/ProjetECO/Cours' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 728 features and 30 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 4.74204 ymin: 44.69597 xmax: 6.359009 ymax: 45.88339
## Geodetic CRS:  WGS 84
sirene_shp <- st_read(
  dsn = here("Cours"),
  layer = "dataseed-sirene-2"
)
## Reading layer `dataseed-sirene-2' from data source 
##   `/Users/loic/Documents/UGA/donnee_entreprise/ProjetECO/Cours' 
##   using driver `ESRI Shapefile'
## replacing null geometries with empty geometries
## Simple feature collection with 342 features and 106 fields (with 5 geometries empty)
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 4.767436 ymin: 44.90643 xmax: 6.121839 ymax: 45.81497
## Geodetic CRS:  WGS 84
iris_data <- read_excel(here("Cours", "iris_isere.xlsx"))

# Zone d'étude: Grenoble + communes nord
sirene_coms <- c(
  "GRENOBLE", "MEYLAN", "LA TRONCHE",
  "CORENC", "SAINT-MARTIN-LE-VINOUX", "SAINT-EGREVE"
)
iris_coms <- c(
  "Grenoble", "Meylan", "La Tronche",
  "Corenc", "Saint-Martin-le-Vinoux", "Saint-Égrève"
)

sirene_grenoble <- sirene_shp %>%
  filter(libellecomm %in% sirene_coms & etatadminis == "Actif")
iris_grenoble <- iris_shp %>%
  filter(com_name %in% iris_coms)

# Préparation des variables IRIS (même logique que test.Rmd)
colnames(iris_grenoble)[20] <- "CODE_IRIS"

iris_data$nb_residence_princ <- iris_data$Res_princ_30_moins_40_m2 +
  iris_data$Res_princ_40_moins_60_m2 +
  iris_data$Res_princ_60_moins_80_m2 +
  iris_data$Res_princ_80_moins_100_m2 +
  iris_data$Res_princ_100_moins_120_m2

vars_to_keep <- c(
  "CODE_IRIS",
  "1er_quartile_euro", "3e_quartile_euro",
  "Actifs_15-64_ans2", "Pop_15-64_ans",
  "Delta_Pop_20-64_ans", "nb_residence_princ"
)

iris_grenoble <- merge(
  iris_grenoble,
  iris_data[, vars_to_keep],
  by = "CODE_IRIS",
  all.x = TRUE
)

iris_grenoble <- iris_grenoble %>%
  rename(
    premier_quartile_revenu = `1er_quartile_euro`,
    troisieme_quartile_revenu = `3e_quartile_euro`,
    Actif_total = `Actifs_15-64_ans2`,
    Pop_15_64 = `Pop_15-64_ans`,
    Delta_pop_20_64ans = `Delta_Pop_20-64_ans`,
    Nb_residence_principale = nb_residence_princ
  )

iris_grenoble <- iris_grenoble %>%
  mutate(
    across(
      c(
        premier_quartile_revenu,
        troisieme_quartile_revenu,
        Actif_total,
        Pop_15_64,
        Delta_pop_20_64ans,
        Nb_residence_principale
      ),
      ~ na_if(., 0)
    )
  )

st_crs(sirene_grenoble)
## Coordinate Reference System:
##   User input: WGS 84 
##   wkt:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["latitude",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["longitude",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     ID["EPSG",4326]]
# Buffer de 500 m autour des supermarchés
sirene_grenoble_l93 <- st_transform(sirene_grenoble, 2154)
buffer <- st_buffer(sirene_grenoble_l93, 500)
buffer <- st_transform(buffer, st_crs(iris_grenoble))

tmap_mode("view")
tm_shape(iris_grenoble) +
  tm_text(text = "iris_name", size = 0.8) +
  tm_polygons(
    "Pop_15_64",
    alpha = 0.5,
    style = "pretty",
    id = "iris_name",
    title = "Population (15-64 ans)"
  ) +
  tm_shape(buffer) +
  tm_borders("chartreuse3", lwd = 3) +
  tm_shape(sirene_grenoble) +
  tm_dots(
    "enseigne1et",
    legend.show = TRUE,
    title = "Enseigne",
    id = "enseigne1et",
    popup.vars = "trancheeffe.2",
    alpha = 0.3
  ) +
  tm_layout(legend.outside = TRUE)
tm_shape(iris_grenoble) +
  tm_text(text = "iris_name", size = 0.8) +
  tm_polygons(
    "Pop_15_64",
    alpha = 0.5,
    style = "pretty",
    id = "iris_name",
    title = "Population (15-64 ans)"
  ) +
  tm_shape(buffer) +
  tm_polygons("enseigne1et", legend.show = FALSE, id = "enseigne1et", alpha = 0.5) +
  tm_shape(sirene_grenoble) +
  tm_dots(
    "enseigne1et",
    legend.show = TRUE,
    title = "Enseigne",
    id = "enseigne1et",
    popup.vars = "trancheeffe.2",
    alpha = 0.8
  ) +
  tm_layout(legend.outside = TRUE)
## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_polygons()`: instead of `style = "pretty"`, use fill.scale =
## `tm_scale_intervals()`.
## ℹ Migrate the argument(s) 'style' to 'tm_scale_intervals(<HERE>)'
## [v3->v4] `tm_polygons()`: use `fill_alpha` instead of `alpha`.
## [v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
## [v3->v4] `tm_polygons()`: use `fill_alpha` instead of `alpha`.
## [v3->v4] `tm_polygons()`: use `fill.legend = tm_legend_hide()` instead of
## `legend.show = FALSE`.
## [v3->v4] `tm_dots()`: use `fill_alpha` instead of `alpha`.
## [tm_dots()] Arguments `legend.show` and `title` unknown.

2. Obtenir des isochrones

2.1 Mise en place du calculateur d’itinéraire

  • Nous allons utiliser OpenTripPlanner, un logiciel codé en Java et utilisé notamment dans le calculateur d’itinéraire du TAG. Pour fonctionner, OTP a besoin de 3 choses:

Un dossier & un sous-dossier specifiques à créer dans votre R working directory > OTP > graphs > default

Des données sur le réseau de transport GTFS (General Transit Feed Specification) à déposer dans OTP > graphs > default

Un fond de carte openstreetmap à déposer dans OTP > graphs > default

  • données GTFS
gtfs_img_candidates <- c(
  here("Cours", "gtfs.png"),
  here("Cours", "Capture.PNG"),
  here("Cours", "Capture.png")
)
gtfs_img <- gtfs_img_candidates[file.exists(gtfs_img_candidates)][1]
if (!is.na(gtfs_img)) {
  knitr::include_graphics(gtfs_img)
} else {
  cat("Ajoutez une capture GTFS dans Cours/gtfs.png pour l'affichage.")
}

Vous pouvez par exemple trouver les données GTFS sur le site M mobilité du TAG: data.mobilites-m.fr
Une fois téléchargé le fichier SEM-GTFS.zip renommez le en gtfs.zip et placez le dans OTP > graphs > default
  • fond de carte OSM
Il y a aussi beaucoup de sources, par exemple :
https://download.geofabrik.de/europe/france/rhone-alpes.html
Placez le fichier dans OTP > graphs > default


  • Mise en place d’OTP
path_data <- NULL
for (candidate in c(here("Cours", "OTP"), here("Cours", "data"))) {
  if (
    dir.exists(file.path(candidate, "graphs", "default")) ||
    file.exists(file.path(candidate, "Graph.obj"))
  ) {
    path_data <- candidate
    break
  }
}
if (is.null(path_data)) {
  path_data <- here("Cours", "OTP")
  dir.create(path_data, recursive = TRUE, showWarnings = FALSE)
}
path_otp <- otp_dl_jar()

# Vous n'avez besoin de lancer la construction du graphe qu'une fois
# log <- otp_build_graph(otp = path_otp, dir = path_data, memory = 15000)

log1 <- otp_setup(otp = path_otp, dir = path_data)
otpcon <- otp_connect()
otp_img_candidates <- c(
  here("Cours", "OTP_capture.png"),
  here("Cours", "OTP_capture.PNG")
)
otp_img <- otp_img_candidates[file.exists(otp_img_candidates)][1]
if (!is.na(otp_img)) {
  knitr::include_graphics(otp_img)
} else {
  cat("Ajoutez une capture OTP dans Cours/OTP_capture.png pour l'affichage.")
}

2.2 Estimer les isochrones (otp_isochrone())

iso <- otp_isochrone(
  otpcon,
  fromPlace = sirene_grenoble,
  fromID = sirene_grenoble$siret,
  mode = c("WALK", "TRANSIT"),
  maxWalkDistance = 2000,
  date_time = as.POSIXct(strptime("2022-02-22 08:35", "%Y-%m-%d %H:%M")),
  cutoffSec = c(5, 7, 9, 11) * 60
)
iso$minutes <- iso$time / 60

# Vérification
colnames(iso)[3] <- "siret"
iso <- merge(
  iso,
  st_drop_geometry(sirene_grenoble[, c("siret", "enseigne1et")]),
  by = "siret",
  all.y = TRUE
) %>%
  unique()

# Option de correction ponctuelle (si un point ne retourne pas d'isochrone)
if (nrow(sirene_grenoble) >= 16) {
  st_geometry(sirene_grenoble)[16] <- st_sfc(
    st_point(c(5.728253270386639, 45.17659752397455)),
    crs = st_crs(sirene_grenoble)
  )
}

# On recommence
iso <- otp_isochrone(
  otpcon,
  fromPlace = sirene_grenoble,
  fromID = sirene_grenoble$siret,
  mode = c("WALK", "TRANSIT"),
  maxWalkDistance = 1500,
  date_time = as.POSIXct(strptime("2022-02-22 08:35", "%Y-%m-%d %H:%M")),
  cutoffSec = c(5, 7, 9, 11) * 60
)
colnames(iso)[3] <- "siret"

iso <- merge(
  iso,
  st_drop_geometry(sirene_grenoble[, c("siret", "enseigne1et")]),
  by = "siret",
  all.y = TRUE
) %>%
  unique()
iso$minutes <- iso$time / 60

3. Cartes finales et indices

3.1 Carte intéractive avec isochrones

tmap_mode("view")
tm_shape(iris_grenoble) +
  tm_polygons(
    "Pop_15_64",
    alpha = 0.5,
    style = "pretty",
    id = "iris_name",
    title = "Population (15-64 ans)"
  ) +
  tm_shape(sirene_grenoble) +
  tm_dots(
    "enseigne1et",
    legend.show = TRUE,
    title = "Enseigne",
    size = 0.1,
    id = "enseigne1et",
    popup.vars = "trancheeffe.2"
  ) +
  tm_layout(legend.outside = TRUE) +
  tm_shape(iso) +
  tm_fill(
    "minutes",
    breaks = c(0, 5.01, 7.01, 9.01, 11.01),
    title = "Isochrones (minutes)",
    style = "fixed",
    labels = c("0 to 5", "5 to 7", "7 to 9", "9 to 11"),
    palette = "-BuPu",
    id = "minutes",
    alpha = 0.3
  ) +
  tm_borders()

3.2 Analyse de la concurrence

  • Calculez le nombre de supermarchés accessibles par IRIS (st_intersects())
iso420 <- iso %>% filter(time == 420)

iris_grenoble$nb_supermarche_per_iris <- rowSums(
  ifelse(st_intersects(x = iris_grenoble, y = iso420, sparse = FALSE), 1, 0)
)

tm_shape(iris_grenoble) +
  tm_polygons(
    "nb_supermarche_per_iris",
    style = "pretty",
    id = "iris_name",
    labels = c("0", "1", "2", "3", "4", "5"),
    title = "Nb de Supermarchés accessible en 7min"
  ) +
  tm_text(text = "iris_name", size = 0.8)
  • Bonus pour la 3ème session, on calcule le nombre de supermarchés de chaque enseigne accessibles par IRIS
data_temp <- matrix(0, nrow = nrow(iris_grenoble), ncol = nrow(iso420))
for (j in 1:nrow(iso420)) {
  for (i in 1:nrow(iris_grenoble)) {
    data_temp[i, j] <- as.character(st_intersects(iris_grenoble, iso420, sparse = FALSE)[i, j])
    data_temp[i, j] <- ifelse(data_temp[i, j] == "TRUE", iso420$enseigne1et[j], 0)
  }
}

data_temp2 <- matrix(0, nrow = nrow(data_temp), ncol = length(unique(iso420$enseigne1et)))
for (j in 1:length(unique(iso420$enseigne1et))) {
  data_temp2[, j] <- apply(
    data_temp,
    1,
    function(x) length(which(x == unique(iso420$enseigne1et)[j]))
  )
}

data_temp2 <- as.data.frame(data_temp2)
colnames(data_temp2) <- unique(iso420$enseigne1et)
write.csv(data_temp2, here("Cours", "nb_pdv_iris.csv"), row.names = FALSE)
  • Calculez la surface cumulée couverte par un supermarché pour chaque iris (st_intersection(), st_area())
intersect_pct <- st_intersection(iris_grenoble, iso420) %>%
  mutate(intersect_area = st_area(.)) %>%
  dplyr::select(iris_name, intersect_area) %>%
  group_by(iris_name) %>%
  mutate(intersect_area = sum(intersect_area)) %>%
  st_drop_geometry() %>%
  unique()
  • Puis calculez la surface de chaque IRIS, le taux de couverture par iris et faîtes une carte
iris_grenoble <- mutate(iris_grenoble, iris_area = st_area(iris_grenoble))

iris_grenoble <- merge(iris_grenoble, intersect_pct, by = "iris_name", all.x = TRUE)

iris_grenoble <- iris_grenoble %>%
  mutate(intersect_area = ifelse(is.na(intersect_area), 0, intersect_area)) %>%
  mutate(couverture = as.numeric(intersect_area / iris_area) * 100)

tm_shape(iris_grenoble) +
  tm_polygons(
    col = "couverture",
    alpha = 0.5,
    style = "cont",
    id = "iris_name",
    title = "Couverture (%)"
  ) +
  tm_text(text = "iris_name", size = 0.8) +
  tm_shape(sirene_grenoble) +
  tm_dots(
    "enseigne1et",
    legend.show = TRUE,
    title = "Enseigne",
    size = 0.1,
    id = "enseigne1et",
    popup.vars = "trancheeffe.2"
  ) +
  tm_layout(legend.outside = TRUE)
  • Enregistrer iris_grenoble pour la prochaine séance
st_write(
  iris_grenoble,
  here("Cours", "iris_grenoble2.gpkg"),
  delete_dsn = TRUE
)

Devoirs pour la prochaine séance

Préparez plusieurs cartes avec des isochrones différentes et d’autres variables au niveau des IRIS (carte 3.1) afin de proposer une localisation potentielle d’un nouveau supermarché à Grenoble. Proposez aussi une autre carte (3.2) construite à partir d’une isochrone différente (>7min). Ce travail sera évalué.